//see explanation in 
//int main()

#include <fstream>

#include "Math/Functor.h"
#include "Math/Factory.h"
#include "Math/Minimizer.h"

#include "TASImage.h"
#include "TCanvas.h"
#include "TColor.h"
#include "TChain.h"
#include "TF1.h"
#include "TFile.h"
#include "TGaxis.h"
#include "TGraph.h"
#include "TGraphAsymmErrors.h"
#include "TGraphErrors.h"
#include "TH1D.h"
#include "TH2D.h"
#include "TLine.h"
#include "TMath.h"
#include "TMatrixD.h"
#include "TMinuit.h"
#include "TLatex.h"
#include "TLegend.h"
#include "TPaletteAxis.h"
#include "TRandom3.h"
#include "TROOT.h"
#include "TStopwatch.h"
#include "TStyle.h"
#include "TTree.h"
#include "TVector3.h"

#define EPSILON 1e-12

typedef Double_t (*FFunc)(const Double_t *xx, const Double_t *par);

enum opt{
  kTMINUIT=0,
  kMINUIT=1,
  kMINUIT2=2
};

Int_t fgFitPrintLevel = 1;
Int_t fgFitNStep=0;
Int_t fgFitNData=0;
Double_t * fgFitX = 0x0;
Double_t * fgFitY = 0x0;
Double_t * fgFitEx = 0x0;
Double_t * fgFitEy = 0x0;
FFunc fgFitFunc = 0x0;
Double_t fTol = 1;
Int_t fMaxIter = 1e8;
Int_t fMinOpt = kTMINUIT;

Int_t fAlephOpt = -999;

Double_t GetFCNChisquare(const Double_t pars[])
{
  Double_t f = 0;

  for(Int_t idata=0; idata<fgFitNData; idata++){
    const Double_t evaly = fgFitFunc(&(fgFitX[idata]), pars);
    const Double_t dy = fgFitY[idata] - evaly;
    
    Double_t dfdx=0;
    //calculated df/dx if necessary
    if( fgFitEx[idata]!=0){
      const Double_t derix0=fgFitX[idata]-fgFitEx[idata];
      const Double_t derix1=fgFitX[idata]+fgFitEx[idata];
      const Double_t deriy0=fgFitFunc(&derix0, pars);
      const Double_t deriy1=fgFitFunc(&derix1, pars);
      dfdx = (deriy1-deriy0)/(derix1-derix0);
    }
    f += dy*dy/(fgFitEy[idata]*fgFitEy[idata] + dfdx*fgFitEx[idata]*dfdx*fgFitEx[idata]);
  }

  fgFitNStep++;

  return f;
}

void MinFCNChisquare(Int_t &npar, Double_t *gin, Double_t &f, Double_t *pars, Int_t iflag)
{
  f = GetFCNChisquare(pars);
}

Int_t MinimizerFitKernel(const std::string minName, const Bool_t kml,const Int_t npar, Double_t pars[], Double_t errs[], Double_t chi[], const Bool_t *pfix, Double_t *cov, const Double_t *lows, const Double_t *highs)
{
  TStopwatch stwa;
  if(fgFitPrintLevel>=0){
    stwa.Start(kTRUE);
  }

  std::string algoName;
  algoName="Migrad";

  ROOT::Math::Minimizer *minimizer = ROOT::Math::Factory::CreateMinimizer(minName, algoName); 
  
  ROOT::Math::Functor func(/*kml ? &GetFCNMaximumLikelihood :*/ &GetFCNChisquare, npar);
  minimizer->SetFunction(func); 
  
  minimizer->SetPrintLevel(fgFitPrintLevel); 
  
  const Double_t errdef = kml? 0.5 : 1;
  minimizer->SetErrorDef(errdef);

  minimizer->SetStrategy(1);

  minimizer->SetMaxFunctionCalls(fMaxIter);
  minimizer->SetMaxIterations(fMaxIter);

  const Double_t tol = (minName=="Minuit2"?0.5:1)*errdef*fTol;//0.5 due to internal definition of edm in Minuit2
  minimizer->SetTolerance(tol);
  
  Int_t nfree=0;
  for(Int_t ipar=0;ipar<npar;ipar++){
    const Double_t start=pars[ipar];
    const Double_t ss=1e-3; //TMath::Abs(pars[ipar])>EPSILON ? TMath::Abs(pars[ipar])*1e-1 : 1e-3;
    const Double_t low= (lows && lows[ipar]>-1e9)?lows[ipar]:0;
    const Double_t hi= (highs && highs[ipar]>-1e9)?highs[ipar]:0;
 
    if(pfix && pfix[ipar]){
      minimizer->SetFixedVariable(ipar, Form("p%d",ipar), start);
    }
    else{
      if(low<hi){
        minimizer->SetLimitedVariable(ipar, Form("p%d",ipar), start, ss, low, hi);
      }
      else{
        minimizer->SetVariable(ipar, Form("p%d",ipar), start, ss);
      }
      nfree++;
    }
  }
  
  //Minuit2-Migrad
  //http://root.cern.ch/root/html/ROOT__Minuit2__Minuit2Minimizer.html
  //status = 1    : Covariance was made pos defined
  //status = 2    : Hesse is invalid
  //status = 3    : Edm is above max
  //status = 4    : Reached call limit
  //status = 5    : Any other failure

  minimizer->Minimize(); 
  Int_t migradflag = minimizer->Status();
  
  //get par right after MIGRAD before HESSE modifies them
  Double_t migp[10];
  memcpy(migp, minimizer->X(), sizeof(Double_t)*npar);

  minimizer->Hesse();
  const Int_t hesseflag = minimizer->Status()-migradflag;

  //only correct migradflag after hesseflag set
  if(minName=="Minuit2" && migradflag==1){
    //pos def not required here
    migradflag=0;
  }

  const Int_t kfail = migradflag+hesseflag;
  if(kfail){
    printf("FitKernel fitting error !! %d %d!!!\n", migradflag, hesseflag);
  }

  if(errs){
    memcpy(errs, minimizer->Errors(), sizeof(Double_t)*npar);
  }

  memcpy(pars, migp, sizeof(Double_t)*npar);
  if(chi){
    chi[0] = GetFCNChisquare(pars);
    chi[1] = -nfree;//not working: minimizer->NFree();//need to add nx outside FitKernel
  }

  if(cov){
    minimizer->GetCovMatrix(cov);
  }

  delete minimizer;

  if(fgFitPrintLevel>=0){
    stwa.Stop(); 
    stwa.Print("m"); //"m": milli sec; "u": micro sec
  } 

 //test->
  for(Int_t ii=0; ii<npar; ii++){
    printf("test0 par %d %e\n", ii, pars[ii]);
    printf("test1 err %d %e\n", ii, errs[ii]);
  }
  //<-test

 return kfail;
}

Int_t TMinuitFitKernel(const Bool_t kml,const Int_t npar, Double_t pars[], Double_t errs[], Double_t chi[], const Bool_t *pfix, Double_t *cov, const Double_t *lows, const Double_t *highs)
{
  //
  //fitkernel for both ML and chi2 fit
  //

  TStopwatch stwa;
  if(fgFitPrintLevel>=0){
    stwa.Start(kTRUE);
  }
  //====>

  TMinuit * mnt = new TMinuit(npar);
  mnt->SetFCN(/* kml ? MinFCNMaximumLikelihood : */ MinFCNChisquare);
  mnt->SetPrintLevel(fgFitPrintLevel);
    
  //http://wwwasdoc.web.cern.ch/wwwasdoc/WWW/minuit/minmain/chapter2_8.html#SECTION0081000000000000000
  //SET ERRordef<up>
  //Sets the value of UP (default value= 1.), defining parameter errors. Minuit defines parameter errors as the change in parameter value required to change the function value by UP. Normally, for chisquared fits UP=1, and for negative log likelihood, UP=0.5. 
  const Double_t errdef = kml? 0.5 : 1;
  mnt->SetErrorDef(errdef);
    
  Int_t nfree = 0;
  for(Int_t ipar=0;ipar<npar;ipar++){
    const Double_t start=pars[ipar];
    //1e-3*par seems too small
    const Double_t ss= TMath::Abs(pars[ipar])>EPSILON ? TMath::Abs(pars[ipar])*1e-1 : 1e-3;//important if limits are set
    const Double_t ll=(lows && lows[ipar]>-1e9)?lows[ipar]:0;
    const Double_t hh=(highs && highs[ipar]>-1e9)?highs[ipar]:0;
    
    mnt->DefineParameter(ipar, Form("p%d",ipar), start, ss, ll,hh);
    nfree++;
    if(pfix && pfix[ipar]){
      mnt->FixParameter(ipar);
      nfree--;
    }
  }

  //http://wwwasdoc.web.cern.ch/wwwasdoc/WWW/minuit/minmain/chapter2_8.html#SECTION0081000000000000000
  //MIGrad[maxcalls] [tolerance]
  //Causes minimization of the function by the method of Migrad, the most efficient and complete single method, recommended for general functions (see also MINImize). The minimization produces as a by-product the error matrix of the parameters, which is usually reliable unless warning messages are produced. The optional argument [maxcalls] specifies the (approximate) maximum number of function calls after which the calculation will be stopped even if it has not yet converged. The optional argument [tolerance] specifies required tolerance on the function value at the minimum. The default tolerance is 0.1, and the minimization will stop when the estimated vertical distance to the minimum (EDM) is less than 0.001*[tolerance]*UP (see [SET ERRordef]SET ERR). 

  Double_t arg[2]={fMaxIter, errdef*fTol};
  Int_t migradflag = -999;
  mnt->mnexcm("MIGRAD", arg, 2, migradflag);

  //get par right after MIGRAD before HESSE modifies them, initial pars not changed
  Double_t migp[10];
  Double_t dummy=-999;
  for(Int_t ipar=0;ipar<npar;ipar++){
    mnt->GetParameter(ipar,migp[ipar], dummy);
  }

  Int_t hesseflag = -999;
  mnt->mnexcm("HESSE", arg ,1, hesseflag);

  const Int_t kfail = migradflag+hesseflag;
  if(kfail){
    printf("FitKernel fitting error !! %d %d\n", migradflag, hesseflag);
  }
  
  for(Int_t ipar=0;ipar<npar;ipar++){
    if(errs) {
      mnt->GetParameter(ipar, dummy,errs[ipar]);
    }
    
    pars[ipar]=migp[ipar];
  }
  if(chi){
    chi[0] = GetFCNChisquare(pars);
    if(nfree!=mnt->GetNumFreePars()){
      printf("TMinuitFitKernel nfree!=mnt->GetNumFreePars() %d %d\n", nfree, mnt->GetNumFreePars());exit(1);
    }
    chi[1] = -nfree;//need to add nx outside FitKernel
  }

  if(cov){
    mnt->mnemat(cov,npar);
  }

  delete mnt;

  //<===
 
  if(fgFitPrintLevel>=0){
    stwa.Stop(); 
    stwa.Print("m"); //"m": milli sec; "u": micro sec
  } 
 
  //test->
  for(Int_t ii=0; ii<npar; ii++){
    printf("test0 par %d %e\n", ii, pars[ii]);
    printf("test1 err %d %e\n", ii, errs[ii]);
  }
  //<-test

  return kfail;
}

Int_t FitKernel(const Bool_t kml,const Int_t npar, Double_t pars[], Double_t errs[], Double_t chi[], const Bool_t *pfix, Double_t *cov, const Double_t *lows, const Double_t *highs)
{
  static Int_t counter = 0;
  if(counter==0){
    printf("\n\n**************** FitKernel fMinOpt %d **************** \n", fMinOpt);
  }

  counter++;

  if(fMinOpt==kTMINUIT){
    return TMinuitFitKernel(kml, npar, pars, errs, chi, pfix, cov, lows, highs);
  }
  else{
    std::string minName;
    if(fMinOpt==kMINUIT){
      minName="Minuit";//this is exactly the same as TMinuitFitKernel
    }
    else if(fMinOpt==kMINUIT2){
      minName="Minuit2";
    }
    else{
      printf("FitKernel fMinOpt error!! %d\n", fMinOpt); exit(1);
    }
    return MinimizerFitKernel(minName, kml, npar, pars, errs, chi, pfix, cov, lows, highs);
  }
}

Int_t ChisquareFit(const Int_t nx, const Double_t xdata[], const Double_t ydata[], const Double_t *xerr, const Double_t yerr[], FFunc ffunc, const Int_t npar, Double_t pars[], Double_t *errs=0x0, Double_t *chi=0x0, const Bool_t *pfix=0x0, Double_t *cov=0x0, const Double_t *lows=0x0, const Double_t *highs=0x0)
{
  if(nx<=0){
    printf("ChisquareFit error nx<=0 %d\n", nx); exit(1);
  }

  fgFitNStep=0;
  fgFitNData = nx;
  fgFitX = new Double_t[nx];
  fgFitY = new Double_t[nx];
  fgFitEx = new Double_t[nx];
  fgFitEy = new Double_t[nx];
  for(Int_t ix=0; ix<fgFitNData; ix++){
    fgFitX[ix]=xdata[ix];
    fgFitY[ix]=ydata[ix];
    fgFitEx[ix]=xerr?xerr[ix]:0;
    fgFitEy[ix]=yerr[ix];
    if(fgFitPrintLevel>=0){
      printf("ChisquareFit %d: %e %e %e %e\n", ix, fgFitX[ix], fgFitY[ix], fgFitEx[ix], fgFitEy[ix]);
    }
  }
  fgFitFunc = ffunc;

  const Int_t kfail = FitKernel(kFALSE, npar, pars, errs, chi, pfix, cov, lows, highs);

  if(chi){
    chi[1] += nx;
    printf("\tChisquareFit Chi2/NDOF %.1f/%.0f = %.1f\n\n", chi[0], chi[1], chi[1]? chi[0]/chi[1] : -999);
  }

  delete fgFitX;
  delete fgFitY;
  delete fgFitEx;
  delete fgFitEy;

  return kfail;
}

Double_t myBetheBloch(const Double_t *xx, const Double_t *par)
{
  //ALEPH has numerical difficulty in ln(p2+1/pow(bg, p4))
  //use flat dedx for electron. deviation from this at high bg is not inportant because anyway there the statistics is very limited
  //npar = 5;

  const Double_t bg = xx[0];
  const Double_t beta = bg/TMath::Sqrt(1.+ bg*bg);

  const Double_t p0 = par[0];//overall prefactor
  const Double_t p1 = par[1];//beta exponent ~ 2
  const Double_t p2 = par[2];
  const Double_t p3 = par[3];
  const Double_t p4 = par[4];//plateau ~ 80

  if(bg<800){
    return p0 *( 1 + p2*TMath::Log(bg) + p3*TMath::Log(beta) )/TMath::Power(beta, p1);
  }
  else{
    return p4;
  }
}

Double_t Lund(const Double_t * xx,  const Double_t * rawpar)
{
  Double_t par[100];
  memcpy(&(par[1]), rawpar, sizeof(Double_t)*5);

  const Double_t bg = xx[0];

  // bg is beta*gamma
  const Double_t beta2 = bg*bg / (1.0+bg*bg);
  const Double_t a = par[1];
  const Double_t b = par[2];
  const Double_t c = par[3];
  const Double_t e = par[4];
  const Double_t f = par[5];
  const Double_t d = TMath::Exp(c*(a-f)/b);
  
  const Double_t powbg = TMath::Power(1.0+bg, c);
  
  const Double_t value = a/TMath::Power(beta2,e) +
    b/c*TMath::Log(powbg/(1.0 + d*powbg));
  return value;
}

Double_t ALEPH(const Double_t * xx,  const Double_t * par)
{
  const Double_t bg = xx[0];
  const Double_t beta = bg/TMath::Sqrt(1.+ bg*bg);

  const Double_t p0 = TMath::Abs(par[0]);
  const Double_t p1 = TMath::Abs(par[1]);
  const Double_t p2 = TMath::Abs(par[2]);
  const Double_t p3 = TMath::Abs(par[3]);
  const Double_t p4 = TMath::Abs(par[4]);

  //return p0*((p1-TMath::Log(p2+1/TMath::Power(bg,p4)))/TMath::Power(beta,p3) - 1 );
  return p0*( (p1- TMath::Log(p2)- TMath::Log(1+1/(TMath::Power(bg,p4)*p2))  )/TMath::Power(beta,p3) - 1 );

  //ALEPH has numerical difficulty in ln(p2+1/pow(bg, p4))
  //use flat dedx for electron. deviation from this at high bg is not inportant because anyway there the statistics is very limited
  /*
  //npar = 5;

  const Double_t bg = xx[0];
  const Double_t beta = bg/TMath::Sqrt(1.+ bg*bg);

  const Double_t p0 = par[0];//overall prefactor
  const Double_t p1 = par[1];//beta exponent ~ 2
  const Double_t p2 = par[2];
  const Double_t p3 = par[3];
  const Double_t p4 = par[4];//plateau ~ 80

  const Double_t turingp = 800;// ~ 100 GeV pion ~ 0.4 GeV electron
  if(bg<turingp){
    return p0 *( 1 + p2*TMath::Log(bg) + p3*TMath::Log(beta) )/TMath::Power(beta, p1);
  }
  else{
    return p4;
  }
  */
  
  //p1 shift not needed!!
  //npar = 6;

  /*
  const Double_t bg = xx[0];
  const Double_t beta = bg/TMath::Sqrt(1.+ bg*bg);
 const Double_t p0 = par[0];//overall prefactor
  const Double_t p1 = par[1];//beta exponent ~ 2
  const Double_t p2 = par[2];
  const Double_t p3 = par[3];
  const Double_t p4 = par[4];//plateau ~ 80

  return p0/TMath::Power(beta,2*p4) + p1/p2*TMath::Log(TMath::Power(1+bg,p2)/(1+p3*TMath::Power(1+bg, p2)));
  */

  /*
  const Double_t p0 = par[0];//overall prefactor
  const Double_t p1 = par[1];//shift ~ -beta^p3/beta^p3 * normalization
  const Double_t p2 = par[2];//beta exponent ~ 2
  const Double_t p3 = par[3];
  const Double_t p4 = par[4];
  const Double_t p5 = par[5];//plateau ~ 80

  const Double_t turingp = 800;// ~ 100 GeV pion ~ 0.4 GeV electron
  if(bg<turingp){
    return p0*((p1-TMath::Log(p2+1/TMath::Power(bg,p4)))/TMath::Power(beta,p3)-1);
  }
  else{
    return p5;
  }
  */
  //
  //ALEPH parametrization for dEdx
  //npar = 5
  //

  /*
  const Double_t bg = xx[0];
  const Double_t beta = bg/TMath::Sqrt(1.+ bg*bg);
  const Double_t gamma = bg/beta;

  const Double_t p0 = TMath::Abs(par[0]);
  const Double_t p1 = TMath::Abs(par[1]);
  const Double_t p2 = TMath::Abs(par[2]);
  const Double_t p3 = TMath::Abs(par[3]);
  const Double_t p4 = TMath::Abs(par[4]);
  const Double_t p5 = TMath::Abs(par[5]);

  const Double_t aa = TMath::Power(beta, p3);

  const Double_t turingp = 800;
  if(bg<turingp){
    if(fAlephOpt==0){  
      return p0*(p1-aa-TMath::Log(p2+TMath::Power(bg,-p4)))/aa;
    }
    else if(fAlephOpt==1){
      return p1*( 1 + p2*TMath::Log(beta) + p4*TMath::Log(gamma) )/aa -p0;
    }
    else if(fAlephOpt==3){
      return p1*( 1 + p2*TMath::Log(bg) + p4*TMath::Log(beta) )/aa -p0;
    }
    else if(fAlephOpt==2){
      return p0*( (p1+(1+p2)*TMath::Log(bg)-(2+p3)*TMath::Log(beta))/TMath::Power(beta, 2+p4) - 1 );
    }
    else{
      printf("wrong aleph opt!! %d\n", fAlephOpt);exit(1);
    }
  }
  else{
    return p5;
  }
  */
  //return p0/(p1-1)*(p1-aa-TMath::Log(1+ 1/(p2*TMath::Power(bg,p4)) ))/aa;
  //return p0/(p0-1)*(p1-aa+TMath::Log(1+1/(p2*TMath::Power(bg,p4))))/aa;

  /*
  Double_t bb = 0;
  //0 = 20, checked!
  //1 = 21, not totally
  if(fAlephOpt==0){      //STATUS=CONVERGED FCN=13199.3 FROM HESSE     STATUS=NOT POSDEF     31 CALLS         210 TOTAL
    const Double_t p2 = TMath::Abs(par[2]);
    //numerically very bad when p2~1e-15, bg~1e3, p4~5.
    bb = TMath::Log( p2 + TMath::Power(bg, -p4) );
  }
  else if(fAlephOpt==1){ //STATUS=CONVERGED FCN=13504   FROM HESSE     STATUS=NOT POSDEF     31 CALLS        1885 TOTAL
    const Double_t p2 = TMath::Abs(par[2]);
    bb = TMath::Log( p2 + TMath::Power(bg+1, -p4) );
  }
  else if(fAlephOpt==10){//STATUS=FAILEDnan FCN=13158.1 FROM HESSE     STATUS=NOT POSDEF     39 CALLS       17411 TOTAL 
    //after redefining p2 (<0) -> ln P2
    const Double_t p2 = par[2];
    const Double_t lbg = TMath::Log(bg);
    bb = p2 + TMath::Log( 1 + TMath::Exp(-p2-p4*lbg) );
  }
  else if(fAlephOpt==11){//STATUS=CONVERGED FCN=13504   FROM HESSE     STATUS=OK             31 CALLS        1703 TOTAL

    //after redefining p2 (<0) -> ln P2
    const Double_t p2 = par[2];
    const Double_t lbg = TMath::Log(bg+1);
    bb = p2 + TMath::Log( 1 + TMath::Exp(-p2-p4*lbg) );
    //absorb p2 in p1
    //bb = -p4*TMath::Log(1+1/(bg+1));//TMath::Log( 1 + p2*TMath::Exp(-p4*lbg) );

    //printf("test3 %e %e %e %e %e --- %e %e %e\n", p0, p1, p2, p3, p4, bg+1, lbg, -p2-p4*lbg);
  }
  else if(fAlephOpt==20){//STATUS=CONVERGED FCN=13199.3 FROM HESSE     STATUS=NOT POSDEF     31 CALLS         210 TOTAL
    const Double_t p2 = TMath::Abs(par[2]);
    const Double_t lbg = TMath::Log(bg);
    bb = TMath::Log(p2) + TMath::Log( 1 + TMath::Exp(-TMath::Log(p2)-p4*lbg) );
  }
  else if(fAlephOpt==21){//STATUS=CONVERGED FCN=13504   FROM HESSE     STATUS=OK             31 CALLS        1835 TOTAL
    const Double_t p2 = TMath::Abs(par[2]);
    const Double_t lbg = TMath::Log(bg+1);
    bb = TMath::Log(p2) + TMath::Log( 1 + TMath::Exp(-TMath::Log(p2)-p4*lbg) );
  }
  else{
    printf("bad fAlephOpt!! %d\n", fAlephOpt); exit(1);
  }
 
  return (p1-aa-bb)*p0/aa;
  */
}

TH1D * GetHfit(const TString hname, FFunc func, const Double_t par[], const Double_t xmin, const Double_t xmax)
{
  TH1D * h1=new TH1D(hname, "", 500, xmin, xmax);

  TAxis *axis=h1->GetXaxis();
  const Int_t bins = axis->GetNbins();

  if (xmin<EPSILON){
    printf("xmin < epsilon! %e\n", xmin); exit(1);
  }

  Double_t *new_bins = new Double_t[bins + 1];

  new_bins[0] = xmin;
  const Double_t factor = pow(xmax/xmin, 1./bins);

  for (int i = 1; i <= bins; i++) {
   new_bins[i] = factor * new_bins[i-1];
  }
  axis->Set(bins, new_bins);
  delete [] new_bins;
  
  for(Int_t ii=1; ii<=h1->GetNbinsX(); ii++){
    const Double_t xx = h1->GetBinCenter(ii);
    const Double_t yy = func(&xx, par);
    printf("GetHfit %s %d %e %e\n", hname.Data(), ii, xx, yy);
    h1->SetBinContent(ii, yy);
  }
  return h1;
}

int main(int argc, char * argv[])
{
  /*
    g++ alephFit.C  -O3 -Wall -Werror -I$ROOTSYS/include  -L$ROOTSYS/lib -lCore -lRIO -lHist -lGraf -lGpad -lMinuit  -lMathCore -o alephFit
  
    r=0; m=0; a=0; ./alephFit alephDataRev$r.root $m $a > see-$r\-$m\-$(echo $a | awk '{printf("%02d",$1)}').log
  */

  for(int ii=0; ii<argc; ii++){
    printf("%d: %s\n", ii, argv[ii]);
  }

  if(argc!=3){
    printf("argc!=3 %d\n", argc); return 1;
  }

  const TString fin(argv[1]);
  TFile::Open(fin);
  const Int_t rev = fin.Contains("Rev0") ? 0 : 1;
  TGraphErrors *gr=(TGraphErrors*)gDirectory->Get(Form("grdataRev%d", rev));
  if(!gr){
    printf("no gr!\n"); gDirectory->ls(); return 1;
  }

  //---
  fMinOpt = atoi(argv[2]);
  //fAlephOpt = atoi(argv[3]);

  Double_t par[10];//={100, 2, 0, 0, 80};

  Double_t err[10], chi[10];

  TCanvas *c1=new TCanvas;
  gr->Draw("ap");
  gPad->SetLogx();

  /*
  FFunc func = Lund;//myBetheBloch;//
  par[0] = 34.0446;
  par[1] = 8.42221;
  par[2] = 4.16724;
  par[3] = 1.29473;
  par[4] = 80.6663;
  */

  FFunc func = ALEPH;
  par[0]=       4.401269;
  par[1]=       9.725370;
  par[2]=       0.000178;
  par[3]=       1.904962;
  par[4]=       1.426576;

  const Int_t npar = 5;

  TH1D *hfit0 = GetHfit("prefit", func, par, 0.1, 4e4);
  hfit0->SetLineColor(kBlue);
  hfit0->Draw("same");

  /*
  Double_t lows[]={0,0,0,0,0,0};
  Double_t highs[]={0,0,0,0,0,0};
  Bool_t pfix[]={0,1,0,0,0,0};
  */
  ChisquareFit(gr->GetN(), gr->GetX(), gr->GetY(), 0x0, gr->GetEY(), func, npar, par,err,chi);//, pfix, 0x0, lows, highs); 

  TH1D *hfit1 = GetHfit("postfit", func, par, 0.1, 4e4);
  hfit1->SetLineColor(kRed);
  hfit1->Draw("same");

  TFile *fout=new TFile(Form("outfitRev%d_%d-%d.root", rev, fMinOpt, fAlephOpt),"recreate");
  gr->Write();
  hfit0->Write();
  hfit1->Write();
  c1->Write();
  fout->Save();
  fout->Close();
  delete fout;

  return 1;
}
